﻿Imports FreeImageAPI
Imports OpenTK
Imports OpenTK.Graphics.OpenGL

Module Helper
    Public Sub LoadTexture(filepath As String, GLTextureID As Integer)
        Dim format As FreeImageAPI.FREE_IMAGE_FORMAT = FreeImageAPI.FreeImage.GetFileType(filepath, 0)
        If format = FreeImageAPI.FREE_IMAGE_FORMAT.FIF_UNKNOWN Then Throw New Exception("未知图像格式")
        If FreeImageAPI.FreeImage.FIFSupportsReading(format) Then
            Dim bmp As FreeImageAPI.FIBITMAP = FreeImageAPI.FreeImage.Load(format, filepath, FreeImageAPI.FREE_IMAGE_LOAD_FLAGS.DEFAULT)
            Dim bpp As Integer = FreeImageAPI.FreeImage.GetBPP(bmp)
            Dim type As FreeImageAPI.FREE_IMAGE_COLOR_TYPE = FreeImageAPI.FreeImage.GetColorType(bmp)
            Dim width As Integer = FreeImageAPI.FreeImage.GetWidth(bmp)
            Dim height As Integer = FreeImageAPI.FreeImage.GetHeight(bmp)
            Dim data As IntPtr = FreeImageAPI.FreeImage.GetBits(bmp)

            GL.BindTexture(TextureTarget.Texture2D, GLTextureID)
            GL.TexParameter(TextureTarget.Texture2D, TextureParameterName.TextureMagFilter, TextureMagFilter.Linear)
            GL.TexParameter(TextureTarget.Texture2D, TextureParameterName.TextureMinFilter, TextureMinFilter.Linear)
            Select Case bpp
                Case 24
                    GL.TexImage2D(TextureTarget.Texture2D, 0, PixelInternalFormat.Rgb, width, height, 0, PixelFormat.Bgr, PixelType.UnsignedByte, data)
                Case 32 And type = FreeImageAPI.FREE_IMAGE_COLOR_TYPE.FIC_RGB
                    GL.TexImage2D(TextureTarget.Texture2D, 0, PixelInternalFormat.Rgb, width, height, 0, PixelFormat.Bgra, PixelType.UnsignedByte, data)
                Case 32 And type = FreeImageAPI.FREE_IMAGE_COLOR_TYPE.FIC_RGBALPHA
                    GL.TexImage2D(TextureTarget.Texture2D, 0, PixelInternalFormat.Rgba, width, height, 0, PixelFormat.Bgra, PixelType.UnsignedByte, data)
                Case Else
                    Throw New Exception("未知图像数据")
            End Select
            GL.BindTexture(TextureTarget.Texture2D, 0)
            FreeImageAPI.FreeImage.Unload(bmp)
        End If
    End Sub

    Function ToVector3(v As Assimp.Vector3D) As OpenTK.Vector3
        Return New OpenTK.Vector3(v.X, v.Y, v.Z)
    End Function

    Function ToVector3(c As Assimp.Color3D) As OpenTK.Vector3
        Return New OpenTK.Vector3(c.R, c.G, c.B)
    End Function

    Function ToVector4(v As Assimp.Color4D) As OpenTK.Vector4
        Return New OpenTK.Vector4(v.R, v.G, v.B, v.A)
    End Function

    Public Function ToMatrix4(ByVal m As Assimp.Matrix4x4) As OpenTK.Matrix4
        Dim a(15) As Single
        a(0) = m.A1 : a(4) = m.A2 : a(8) = m.A3 : a(12) = m.A4
        a(1) = m.B1 : a(5) = m.B2 : a(9) = m.B3 : a(13) = m.B4
        a(2) = m.C1 : a(6) = m.C2 : a(10) = m.C3 : a(14) = m.C4
        a(3) = m.D1 : a(7) = m.D2 : a(11) = m.D3 : a(15) = m.D4
        'Return New OpenTK.Matrix4(a(0), a(4), a(8), a(12),
        '                   a(1), a(5), a(9), a(13),
        '                   a(2), a(6), a(10), a(14),
        '                   a(3), a(7), a(11), a(15))
        REM GL的矩阵
        Return New OpenTK.Matrix4(a(0), a(1), a(2), a(3),
                           a(4), a(5), a(6), a(7),
                           a(8), a(9), a(10), a(11),
                           a(12), a(13), a(14), a(15))
    End Function

    Public Function ToMatrix4x4(m As OpenTK.Matrix4) As Assimp.Matrix4x4
        Dim am As Assimp.Matrix4x4
        am.A1 = m.M11 : am.A2 = m.M21 : am.A3 = m.M31 : am.A4 = m.M41
        am.B1 = m.M12 : am.B2 = m.M22 : am.B3 = m.M32 : am.B4 = m.M42
        am.C1 = m.M13 : am.C2 = m.M23 : am.C3 = m.M33 : am.C4 = m.M43
        am.D1 = m.M14 : am.D2 = m.M24 : am.D3 = m.M34 : am.D4 = m.M44

        'am.A1 = m.M11 : am.A2 = m.M12 : am.A3 = m.M13 : am.A4 = m.M14
        'am.B1 = m.M21 : am.B2 = m.M22 : am.B3 = m.M23 : am.B4 = m.M24
        'am.C1 = m.M31 : am.C2 = m.M32 : am.C3 = m.M33 : am.C4 = m.M34
        'am.D1 = m.M41 : am.D2 = m.M42 : am.D3 = m.M43 : am.D4 = m.M44
        Return am
    End Function

    Structure Complex
        Public r As Double
        Public v As Double
        Public Sub New(r As Double)
            Me.r = r
        End Sub

        Public Sub New(r As Double, v As Double)
            Me.r = r
            Me.v = v
        End Sub

        Public Shared Operator +(ByVal a As Complex, ByVal b As Complex) As Complex
            Return New Complex(a.r + b.r, a.v + b.v)
        End Operator

        Public Shared Operator -(ByVal a As Complex, ByVal b As Complex) As Complex
            Return New Complex(a.r - b.r, a.v - b.v)
        End Operator

        Public Shared Operator *(ByVal a As Complex, ByVal b As Complex) As Complex
            Return New Complex(a.r * b.r - a.v * b.v, a.r * b.v + a.v * b.r)
        End Operator

        Public Shared Operator /(ByVal a As Complex, ByVal b As Complex) As Complex
            Dim l = b.r * b.r + b.v * b.v
            Return New Complex((a.r * b.r + a.v * b.v) / l, (a.v * b.r - a.r * b.v) / l)
        End Operator

        Public Shared Widening Operator CType(ByVal s As Double) As Complex
            Return New Complex(s)
        End Operator
    End Structure

    Public Class SortAbsDesc
        Implements IComparer
        Public Function Compare(x As Object, y As Object) As Integer Implements IComparer.Compare
            Dim a = CType(x, Single)
            Dim b = CType(y, Single)
            Return Math.Abs(b) - Math.Abs(a)
        End Function
    End Class

    Public Function Solve_x1p2(a As Single, b As Single, c As Single) As Complex()
        Dim delta = b * b - 4 * a * c
        If delta > 0 Then
            Dim x1 As Complex = -b / (2 * a) + Math.Sqrt(delta) / (2 * a)
            Dim x2 As Complex = -b / (2 * a) - Math.Sqrt(delta) / (2 * a)
            Return {x1, x2}
        Else
            Dim x1 As Complex = -b / (2 * a) + New Complex(0, Math.Sqrt(delta)) / (2 * a)
            Dim x2 As Complex = -b / (2 * a) - New Complex(0, Math.Sqrt(delta)) / (2 * a)
            Return {x1, x2}
        End If
    End Function

    Public Function Clamp(v As Double, min As Double, max As Double) As Double
        If v < min Then v = min
        If v > max Then v = max
        Return v
    End Function

    Public Function Power(v As Double, e As Double) As Double
        If v > 0 Then
            Return Math.Pow(v, e)
        Else
            Return -Math.Pow(-v, e)
        End If
    End Function

    Public Function Solve_2equation_xyz(a As Double, b As Double, c As Double, i As Double, j As Double, k As Double) As Vector3
        Dim x0 = c * j - b * k
        Dim y0 = a * k - c * i
        Dim z0 = b * i - a * j
        If x0 = 0 Or y0 = 0 Or z0 = 0 Then
            Return Vector3.Zero
        Else
            Return New Vector3(x0 / z0, y0 / z0, 1)
        End If
    End Function

    Public Function Solve_x1p3(a As Double, b As Double, c As Double, d As Double) As Complex()
        Dim omega = (-1 + New Complex(0, Math.Sqrt(3))) * 0.5
        Dim omega2 = (-1 * New Complex(0, Math.Sqrt(3)) - 1) * 0.5
        Dim p As Double = (3 * a * c - b * b) / (3 * a * a)
        Dim q As Double = (27 * a * a * d - 9 * a * b * c + 2 * b * b * b) / (27 * a * a * a)
        Dim div3 = 1.0 / 3.0
        Dim sdelta = Math.Pow(q * 0.5, 2) + Math.Pow(p * div3, 3)
        Dim x1, x2, x3 As Complex
        Dim y1, y2, y3 As Complex
        If (sdelta < 0) Then
            Dim r = Math.Sqrt(-Power(p * div3, 3))
            Dim theta = Math.Acos(Clamp(-q / r / 2, -1, 1)) * div3

            Dim k = 2 * Power(r, div3)
            x1 = k * Math.Cos(theta)
            x2 = k * Math.Cos(theta + 120 / 180 * Math.PI)
            x3 = k * Math.Cos(theta + 240 / 180 * Math.PI)
        Else
            Dim delta = Math.Sqrt(sdelta)

            Dim m, n As Double
            m = Power(-0.5 * q + delta, div3)
            n = Power(-0.5 * q - delta, div3)

            x1 = m + n
            x2 = omega * m + omega2 * n
            x3 = omega2 * m + omega * n
        End If
        y1 = x1 - b / (3 * a)
        y2 = x2 - b / (3 * a)
        y3 = x3 - b / (3 * a)
        Return {y1, y2, y3}
    End Function

    Public Function GetRotationOfCoordinateSystem(vx As Vector3, vy As Vector3, vz As Vector3) As Vector4
        'cos(theta) + (1-cos(theta))xx,   (1-cos(theta))xy - sin(theta)z,   (1-cos(theta)xz+sin(theta)y
        '(1-cos(theta))xy + sin(theta)z   cos(theta) + (1-cos(theta))yy     (1-cos(theta)zy - sin(theta)x
        '(1-cos(theta)xz - sin(theta)y    (1-cos(theta)zy + sin(theta)x     cos(theta) + (1-cos(theta))zz
        '=>
        'cos(theta) + (1-cos(theta))xx = vx.x(A)
        '(1-cos(theta))xy + sin(theta)z = vx.y(B)
        '(1-cos(theta)xz - sin(theta)y = vx.z(C)
        '(1-cos(theta))xy - sin(theta)z = vy.x(D)
        'cos(theta) + (1-cos(theta))yy = vy.y(E)
        '(1-cos(theta)zy + sin(theta)x = vy.z(F)
        '(1-cos(theta)xz+sin(theta)y = vz.x(G)
        '(1-cos(theta)zy - sin(theta)x = vz.y(H)
        'cos(theta) + (1-cos(theta))zz = vz.z(I)
        '旋转轴axis(x,y,z) 为单位向量长度为1，所以x*x+y*y+z*z=1, A+E+I = 1 - 2*cos(theta)
        vx.Normalize() : vy.Normalize() : vz.Normalize()
        Dim theta = Math.Acos((vx.X + vy.Y + vz.Z - 1) / 2)
        Dim T = 1 - Math.Cos(theta)
        Dim NA = (vx.Y + vy.X) / 2, NB = (vz.Y + vy.Z) / 2, NC = (vx.Z + vz.X) / 2
        Dim xSqr = (NA * NC) / (NB * T)
        'Dim ySqr = (NA * NB) / (NC * T)
        'Dim zSqr = (NB * NC) / (NA * T)
        Dim x0 = Math.Sqrt(xSqr)
        Dim y0 = NA / (T * x0)
        Dim z0 = NC / (T * x0)
        'If (y0 * z0 * T - NB > Single.Epsilon) Then
        '    x0 = -x0
        '    y0 = -y0
        '    z0 = -z0
        'End If
        Return New Vector4(x0, y0, z0, theta)
    End Function

    Public Function EigenVector(mtx As Matrix3, value As Single()) As Vector3()
        Dim result As New List(Of Vector3)
        For n = 0 To value.Length - 1
            If value(n) = 0 Then
                Console.WriteLine("lambda is zero")
            End If

            'Dim MatKnown = mtx
            'MatKnown.M11 = MatKnown.M11 - value(n)
            'MatKnown.M22 = MatKnown.M22 - value(n)
            'MatKnown.M33 = MatKnown.M33 - value(n)
            '齐次 有同质同源 均匀的意思。如齐次坐标，不同w表示不同位置，都是同源的，均匀的（乘以常数K依然是其中一部分）。
            '特征向量就是 求齐次方程（项的指数次数相同，等式右边为0）组的非零解。其解是线性均匀的
            '解有无数个，是一个基础解系
            '原矩阵M0在代入λ后，得到M1必定非满秩(结论)，不存在确定解，只会得到基础解系
            Dim a = mtx.M11 - value(n), b = mtx.M12, c = mtx.M13
            Dim d = mtx.M21, e = mtx.M22 - value(n), f = mtx.M23
            Dim g = mtx.M31, h = mtx.M32, i = mtx.M33 - value(n)
            Dim r = Helper.Solve_2equation_xyz(a, b, c, d, e, f)
            If (r = Vector3.Zero) Then
                Dim r1 = Helper.Solve_2equation_xyz(d, e, f, g, h, i)
                If (r1 = Vector3.Zero) Then
                    Dim r2 = Helper.Solve_2equation_xyz(d, e, f, g, h, i)
                    If r2 = Vector3.Zero Then
                        Console.WriteLine("Only find is zero vector")
                        'Throw New Exception("Only find is zero vector")
                    End If
                    r = r2
                End If
                r = r1
            End If
            Dim vec = r.Normalized()
            If r <> Vector3.Zero Then
                Console.WriteLine("EigenValue={0},EigenVector={1}", value(n), vec.ToString())
            End If
            result.Add(vec)
            'ax + by + cz = 0
            'dx + ey + fz = 0
            'gx + hy + iz = 0
            '（1）x = (-by - cz) / a
            ' (ae - bd)*y + (af - cd)*z = 0 y/z = (cd -af) / (ae-bd)
            ' (ahy - bgy + aiz - cgz) = 0

        Next
        Return result.ToArray()
    End Function

    Public Function EigenValue(mtx As Matrix3) As Single()
        Dim lambda(2) As Single
        'ax^3+bx^2+cx+d = 0
        'MatA - lambda * MatIdentity = 0
        'det = a*e*i + b*f*g + c*d*h - c*e*g - i*b*d - f*h*a
        'a(ei-hf) - b(di-gf) +c(dh-eg)
        '另 a=>(a-x), e=>(e-x) i=>(i-x)

        '=(a-x)((e-x)(i-x)-hf) - b(d(i-x)-gf) + c(dh-(e-x)g)
        '=(a-x)(ei+xx-ex-ix-hf) - b(di-dx-gf) + c(dh-eg+gx)
        '=(aei+axx-aex-aix-ahf) - (eix+xxx-exx-ixx-hfx)
        '- (bdi-bdx-bgf) + (cdh-ceg+cgx)
        '= -xxx + axx+exx+ixx + (-aex-aix-eix+hfx+bdx+cgx)
        '+ aei-ahf-bdi+bgf+cdh-ceg
        '= -xxx + (a+e+i)xx + (-ae-ai-ei+hf+bd+cg)x + aei-ahf-bdi+bgf+cdh-ceg
        Dim a = mtx.M11, b = mtx.M12, c = mtx.M13
        Dim d = mtx.M21, e = mtx.M22, f = mtx.M23
        Dim g = mtx.M31, h = mtx.M32, i = mtx.M33
        Dim result = Helper.Solve_x1p3(-1, (a + e + i), (-a * e - a * i - e * i + h * f + b * d + c * g),
        (a * e * i - a * h * f - b * d * i + b * g * f + c * d * h - c * e * g))
        Dim l As New List(Of Single)
        For i = 0 To result.Length - 1
            If result(i).v = 0 Then '只需要实根
                l.Add(result(i).r)
            End If
        Next

        'Dim r3 = result(1) + result(2) '互为共轭复数
        'If r3.v = 0 Then
        '    l.Add(r3.r)
        '    Console.WriteLine("use r1+r2")
        'End If
        Return l.ToArray()
    End Function

    Public Sub OrthogonalSchmidt(ByRef v0 As Vector3, ByRef v1 As Vector3, ByRef v2 As Vector3)
        v0.Normalize()
        v1 = v1 - (Vector3.Dot(v0, v1)) * v0
        v1.Normalize()
        v2 = Vector3.Cross(v0, v1)
    End Sub

    Public Function GetMatrixRank(mtx As Matrix3) As Byte
        Dim rmax = 3
        Return rmax
    End Function
End Module
